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^^' Abstract. Lattice-Boltzmann (LB) simulations are a common tool to numerically 

estimate the permeability of porous media. For valuable results, the porous struc- 

^ ■ ture has to be well resolved resulting in a large computational effort as well as 

(4_| \ high memory demands. In order to estimate the permeability of realistic samples, 

lyj ■ it is of importance to not only implement very efficient codes, but also to choose 

O \ the most appropriate simulation setup to achieve accurate results. With the focus 

c/3 ■ on accuracy and computational effort, we present a comparison between different 

^~>| methods to apply an effective pressure gradient, efficient boundary conditions, as 

"^ ■ well as two LB implementations based on pore-matrix and pore-list data structures. 

(N; 

^ ■ 1 Introduction 

m ■ The lattice-Boltzmann (LB) method is a numerical scheme that is able to simulate the 

. . hydrodynamics of fluids with complex inter facial dynamics and boundaries HHH. Its 

ij \ popularity stems from the broad field of possible application and a fair implemen- 

Q ' tation effort compared to other CFD methods. Unlike schemes that are based on a 

discretization of the Navier-Stokes equations, and therefore represent balance equa- 
tions at the continuum level (macroscopic), the LB method represents the dynamics at 
k><( , mesoscopic level by solving the discretized Boltzmann equation. 

j_j ■ There is an increasing interest in the LB method for simulation of flow in complex 

^ geometries since the end of the 1980's [7J when hydrodynamic simulation methods 

were dominated by finite element schemes that solved the Stokes equation mUl. With 
the advent of more powerful computers it became possible to perform detailed sim- 
ulations of flow in artificially generated geometries Iil0| , tomographic reconstructions 
of sandstone samples [11 -15 J, or fibrous sheets of paper |T6^]. An important property 
estimated using the LB method in those geometries is the permeability 1,17] , and since 
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the porous structure has to be well resolved in order to obtain valuable results, a large 
computational effort as well as high memory demands is required. Therefore, it is 
important to develop very efficient simulation paradigms. 

Different alternative simulation setups have been proposed for permeability esti- 
mation which differ in the computational domain setup, how the fluid is driven, or 
how an effective pressure gradient is being estimated. Further possible differences 
include the choice between single relaxation and multirelaxation time lattice Boltz- 
mann implementations or data structures based on a 3D array containing the whole 
discretized simulation volume including rock matrix and pore space (pore-matrix) in 
contrast to data structures limited to a connected list of pore nodes (pore-list) ll6llT8|. 

The current paper focuses on a detailed comparison of some of these possible im- 
plementation details to accurately estimate the permeability of porous media with the 
LB method. We compare the well known D3Q19 single relaxation (LB-BGK) I.5 .19J and 
multirelaxation time (LB-MRT) ||20]| models together with three alternative setups to 
estimate the permeability utilizing an injection channel of variable length (I-Ch), pres- 
sure boundary conditions (p-BC), or a force density applied over the sample (Vp-S). 
The geometries being investigated are a 3D Poiseuille flow in a square pipe and a 
BCC sphere array. While the first one has the advantage of a minimal discretization 
error, the second one more realistically resembles a natural porous medium. We also 
present a comparison of the efficiency of the LB-codes based on the above mentioned 
pore-matrix and pore-list data structures 1.18.1 . 

2 The Lattice-Boltzmann Method 

The discretized lattice-Boltzmann (LB) equation reads 

ni(x + CiAf,f + AO-ni(x,0 = Af ^^S; Jn^(x,0 - nj'^(x,f)) , (2.1) 

7=1 

where x = {x\,xi,x^^ represents a node. The discretization parameters are Ai and 
b.x, while the discrete velocities c; have the dimension A;t/Af. The variable n, is the 
probable number of particles moving with velocity c,. We use a 3D cubic lattice with 
19 discrete velocities c/, i = 1 . . . 19 known as D3Q19 (see Fig.[T]for a visualization) ||4|. 
The term on the right hand side of Eq. (|2.H) is the linearized collision operator, where 
Sjy is the collision matrix also known as scattering matrix and nj'^(x, f ) is the equilib- 
rium distribution [2J. The macroscopic density p(x, i) and velocity u(x, i) are obtained 
from ni{x,t) as 

JV 
p{x,t) = p'^Y^mix^t), (2.2) 

i=l 
N 

p(x,f)u(x,f) = p°Y^ni{x,t)ci, (2.3) 

i=l 
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Figure 1: D3Q19 cubic lattice with lattice vectors Cj. 

where p° is a reference density. The pressure is given by 

p{x,t) =Cs^p{x,t), 

with 

1 fAx\ 

being the speed of sound IH. The kinematic viscosity v is defined as 



V = 



Cs^At /„ T 



2-^-1 
At 



(2.4) 



(2.5) 



(2.6) 



Within the LB method, the single relaxation time (BGK) or multirelaxation time 
(MRT) methods differ only in the w^ay how the terms Sij and n^'^{x, t) are defined. In 
the case of the LB-BGK scheme IlKlll, the matrix S,y = —Sij/r, where t is a unique 
relaxation time and Sjj is the Kroneker delta. The equilibrium distributions 



P 



<V,0=^c,£, 1 + 



U C; 
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(2.7) 



are obtained from a 2^'^ order Taylor approximation of the Maxwell distribution func- 
tion [|21||, with u* = u* (x, t) and p = p{x, t) calculated in absence of an external force as 
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u*(x, t) = u(x, t) using Eqs. (|2.2|) and (|2.3|l . The numbers a;^; are called lattice weights 
and for the D3Q19 lattice they are 

'1/3, forc/ = 0, 
a;, J 1/18, for I Cj 1 = 1, (2.8) 

_l/36, for|ci| = v/2. 

Generally, they depend on the lattice type, dimensions of space and on the number 
of discrete velocities N. In the article of Qian et ah lH a comprehensive overview on 
different lattices is given. 

In the case of the LB-MRT scheme l|20| the linearized collision operator is rewritten 
as 

Y,SrAnj{x,t)-n'^{x,t)) = -M'^ ■ S ■ {rry{x, t) - m^\x, t)) , (2.9) 

i.e., it is implemented in the space of the hydraulic modes of the problem m(x, t) = 
M • n(x, t) instead of the space of the discrete velocities, where M is the transformation 
matrix. Some modes have a hydrodynamic meaning, but some of them do not and 
are used to improve the numerical stability as described in Il20ll22| . The relaxation 
parameters of the modes are given in the diagonal matrix S. Only two values of the 
diagonal elements of S are used to specify the kinematic viscosity, Eq. I|2.6|) . and the 
bulk viscosity. These are labeled t and t\,^\\^, respectively. In summary, 

S = diag(0,l/Tbuik, 1.4, 0,1.2, 0,1.2, 0,1.2, 1/t, 

1.4, 1/t, 1.4, 1/t, 1/t, 1/t, 1.98, 1.98, 1.98), (2.10) 

is assumed for the MRT method. The other values of the diagonal elements are chosen 
to optimize the performance of the algorithm as described in [20[|22| . The components 
of the equilibrium distribution m'^l(x, t), which represent the density and momentum 
are conserved after collision. The remaining ones are assumed to be functions of the 
conserved moments and for D3Q19 are explicitly given in 1.20 J . 

3 Driving the fluid 

3.1 External force in the LB method 

An external force density of the fluid h{x,t) is implemented by adding a term 

(pi{x,t) = cv,^—-^h{x,t) ■ Ci (3.1) 

to the right hand side of Eq. \2.\) . The equilibrium velocity u*(x, t) is then calculated 
using 

N 

p{x,t)u*{x,t)=p°Y,ni{^'t)^u (3.2) 

!=1 
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but the macroscopic velocity u{x,t) has to be redefined as 

^ At 

p{x, f)u(x, t)=p°Y, M^' i) Ci + -^^i^' 0- (3-3) 

An important point to stress here is that this implementation is just a simplified alter- 
native to add an external acceleration, recovering correctly the hydrodynamic macro- 
scopic fields with less computational cost |23|. In all simulations presented in this 
work the force density acts towards the direction of the X3-axis, b = be^. 

3.2 On-site Pressure Boundary Condition 

Pressure boundary conditions can be implemented as introduced by Zou and He 1241 
|25J. Due to the ideal gas equation of state setting the pressure on a specific node 
is equivalent to setting the density, see Eq. jlAj . The on-site pressure BC are used 
to drive the flow by setting constant densities p™ and p""' (p™ > p°^^) at the inlet 
C(l) and outlet C(X3) planes, i.e. the first and last plane in the direction of the flow, 
respectively. Expanding the calculation of the density p, Eq. (I2.2D , and momentum pu, 
Eq. (|231), leads to 

p 

-L- = m + n2 + "3 + ?^4 + J^5 + "6 + "7 + "8 + "9 + J^IO + '^ll + 
P 

+ni2 + "13 + ni4 + ni5 + nig + ni7 + n^s + nig, (3.4) 

-^Uxj = ni - n2 + "7 + «8 + "9 + mo - flu - ni2 - ni3 - tiu, (3.5) 

p° 

— U;t2 = "3 - n4 + n/ - ns + riu - ni2 + nis + riie - my - nis, (3.6) 

■^Ux3 = ns-ne + m- nw + ni3 - nu + «i5 - "16 + "17 - «18- (3.7) 

Out of the four macroscopic values p and u (uxj, Ux2, and Uxg), only three are inde- 
pendent and can be set, because they have to fulfill the mass balance equation. For 
the calculations performed for this work, the components of the velocity u^j = and 
Ux2 = are being set, resulting in a flow perpendicular to the inlet and outlet plane. 
The variables n^, rig, ni3, ni5, and riiy at the plane C(l), and ne, nio, "14, Wi6/ and nig 
at the plane C(X3) (the first ones and last ones with positive and negative component 
in xo, direction, respectively, see Fig. [l]), are not being updated by the streaming step 
within the LB algorithm but they have to be solved in order to obtain the requested 
density (pressure) and velocity on the node. As it is explained above, since p, u^j, and 
Ux2 are already set, it is not possible to also set the velocity u^,, and its value has to be 
solved for on every node in the inlet plane C(l) and outlet plane C{Xs). Therefore, we 
have a linear system of four equations and six unknowns to be solved. To complete 
the system of equations, bounce-back BC are assumed for the non-equilibrium part of 
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(3.8) 



where the terms N^l and Nxl are the transversal momentum corrections on the X3- 
boundary for the distributions propagating in Xi and X2, respectively ||25| . Eq. (|3.8|) 
then reads 

"9 - "14 = ^U;,3 - N^^f , ni3 - "10 = ^u,3 + N^3, (3.9) 



6-^u..3-N-, „,,-nie = ^„ 



"15 - "18 = 737U.,3 - N^|, fin - file = T7^^X3 + K!- 



The on-site pressure BC presented by Kutay et al. fl6\ uses the same equations pre- 
sented here, but lacks the transversal momentum corrections Nxl and N^^. These 
terms are important when velocity gradients are present ||24|. 

4 Simulation Setup 

The computational domain A" is composed of two subsets: Al represents the solid 
nodes and V represents the fluid nodes, with A^ U T' = A" and Al n P = 0. The 
computational domain A" with the dimensions [Xi x X2 x X3], has three zones as pre- 
sented in Fig.m the inlet I, the sample S, and the outlet O with lengths Xj, Xg, and 
Xq, respectively. In the case of using an injection channel the force density b drives 
the fluid in some interval of X3 within the inlet zone I. Cross-sections perpendicular 
to the direction of flow are denoted by 

C(fl) = {x G A" : X3 = fl}. (4.1) 

For instance, the cross-sections C{Xx + 1) and C(X3 — Xq) represent the first and last 
cross-sections within the sample, respectively. The average density in the cross-section 
C{a) is calculated by 

{P)c{a) = y —■ (4.2) 

xe{C(a)r\V) 
The average mass flux through the whole sample is given by 

Q = Y'LPi^)^^A^)^^^f' (4-3) 
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Figure 2: Computational domain X: sample S and two chambers X and C Cross-sections C are perpen- 
dicular to the direction of the flow. 

where |0(x)ux3 (x) is the momentum component in direction of the flow. The perme- 
ability K as defined by Darcy's law is 



Q 



(4.4) 



where A represents the cross sectional area lllTll . The average pressure gradient in 
flow direction ((V/?)j3) is calculated in different manners according to the way how 
the fluid is driven. This issue is explained in more detail below. 



5 Pore-Matrix and Pore-List 

The LB method is typically implemented representing the pore space and the solid 
nodes using a 3D array including the distribution functions Wj and a Boolean vari- 
able to distinguish between a pore and a matrix node. This 3D array is known as 
pore-matrix. During a simulation, at every time step the loop covering the domain 
includes the fluid and the solid nodes. Therefore, if-statements are necessary to dis- 
tinguish whether the node represents a solid node or a fluid node where the collision 
and streaming steps need to be applied. Furthermore, it has to be determined for ev- 
ery node if BC (periodic and no-slip) need to be applied. The advantage of this data 
structure is its straightforward implementation. However, for the simulation of fluid 
flow in porous media with low porosity the drawbacks are high memory demands 
and inefficient loops through the whole simulation domain. 

An alternative data structure, known as pore-list [18], uses a two-dimensional ar- 
ray characterizing the porous structure. It contains the position (pore-position-list) 
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and connectivity (pore-neighbor-list) of the fluid nodes only. It can be generated from 
the original 3D array before the first time step of the simulation, so that only loops 
through the list of pore nodes not comprising any if-statements for the lattice Boltz- 
mann algorithm are required. The time needed to generate and save the pore-list data 
is comparable to the computational time required for one time step of the LB numer- 
ical scheme implemented by the pore-matrix. This alternative approach is slightly 
more complicated to implement, but allows highly efficient simulations of flows in 
geometries with a low porosity. If the porosity becomes too large, however, the addi- 
tional overhead due to the connection matrix reduces the benefits and at some point 
renders the method less efficient than a standard implementation. 

In Tbl. [l] the amount of memory required for both LB implementations is pre- 
sented. Introducing the notation F as the number of fluid nodes and X'^ = Xi X2 X3 
as the total number of nodes, we can estimate if the pore-list implementation saves 
memory space using the inequality 

(32bit)X3 + 19(64bit)x3 + 19(64bit)X2 > 3(32bit)F + 18(32bit)f + 38(64bit)f . (5.1) 

The terms on the left hand side represent the necessary memory for the pore-matrix 
implementation: pore-matrix, distribution function rij, and its buffer, respectively. The 
terms on the right hand side represent the necessary memory for the pore-list im- 
plementation: pore-position-list, pore-neighbor-list, distribution function n,, and its 
buffer, respectively. Neglecting the term representing the Wj-buffer of the pore-matrix 
implementation, because of being order X^, we can estimate that the pore-list imple- 
mentation saves memory when f /X'^ < 0.40, where F /X^ = (p is the porosity of the 
porous medium. 

6 Computational Domain & Setup 

For the simulations six different lengths of the chambers I and O, Xj and Xq are used 
(see Fig. 12]). We restrict ourselves to Xj = Xq and use the lengths 20, 10, 5, 2, 1 and 0. 
A length of represents simulations without chambers. To drive the flow in order to 
generate an effective pressure gradient, we use three different methods: 

Injection Channel (I-Ch): In the cross-sections C(Xj — 3), C(Xj — 2), and C(Xj — 1) 
a constant force density h {At)^ /{p° Ax) = 10^^ is applied and periodic BC act 

Memory Space 
Pore-Matrix Pore-List 

Geometry: Pore-Matrix, integer [Xi, X2, X3] Pore Position-List, integer [F, 3] 

Pore Neighbor-List, integer [F, 18] 
n,: double-precision [Xi,X2,X3, 19] double-precision [F, 19] 

n; buffer: double-precision [Xi,X2, 3, 19] double-precision [F, 19] 

Table 1: Memory demand of implementations using a pore-matrix or a pore-list data structure. F represents 
the number of fluid nodes. 
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Pressure BC 
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Figure 3: Visualization of the different methods to drive the fluid. The sample is represented as a BCC 
sphere array. I-Ch: the marked zone within T is where the fluid is driven by the force density. p-BC: at the 
marked cross-sections C(l) and C(X3) constant densities are being set. Vp-S: inside the sample (marked 
zone) the fluid is driven by the applied force. 



in all directions. The pressure gradient is estimated using 

2(P)c(X3-Xo) - {p)c{Xi+l) 



((Vp)x3) 



{Xs - l)Ax 



(6.1) 



Pressure BC (p-BC): In the cross-sections C(l) and ^(Xs) constant densities p^ / p° = 
1 + 5 X 10"^ and p°^^ / p° = 1 — 5 x 10^^ are being imposed by using on-site pres- 
sure BC. Periodic BC act in directions Xi and X2. The pressure gradient is calcu- 
lated using Eq. (|6.Hl . 



Sample Force Density (Vp-S): A constant force density h {At)^/{p° Ax) = 10^^ acts 
inside the sample S. In this case it is not necessary to estimate the average pres- 
sure gradient because its value is simply given by b. 

Two different samples are being used, the first one is a 3D Poiseuille flow in a 
square pipe of side length d. The sample dimensions are X2 = X3 = d/ Ax + 2, and 
Xg = 50 + Xj + Xq, where the cross-sectional area in this case is A = d^. The two 
extra nodes of X2 and X3 are used to provide solid walls which define the pipe. The 
theoretical permeability for 3D Poiseuille flow in a square pipe is given by p7| 



,th 




64 



j^D 



n=0 



tanhn2n + l) 
(2n + 1)5 



w 



0.03514^2 



(6.2) 



The second sample is composed of a cube filled with a BCC sphere array. The sample 
geometry is adjusted by varying the radius of the spheres. The dimensions are X2 = 
X3 = 50 and Xg = 50 with Xi = X5 + Xj + X^, giving a cross-sectional area of 
A = 2500 (Ax)^. The theoretical value for the permeability is given by I.28..29.I 



.th 



Gnxh' 



(6.3) 
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where a and h represent the side length of the cube and the inverse of the dimen- 
sionless drag, respectively. As it is presented in 130)1 , h is entirely determined by the 
geometry characteristic of the sphere array and can be represented as a series expan- 
sion of the radius ratio 

X=-^, (6.4) 

'max 

where r represents the radius of the spheres and r^^x = V^ d /4 the maximum radius 
that the spheres can have in the BCC array until they touch each other. The relative 
error calculated by 

e{K) = L_^, (6.5) 

is used to qualify the results of our simulations. 



7 Results 

In order to compare the efficiency and accuracy of different implementations to com- 
pute permeabilities, its value is being estimated by applying the alternatives to drive 
the flow presented above, i.e.. Injection Channel (I-Ch), pressure BC (p-BC), and Sam- 
ple Force Density (V/?-S). It is well know that a slight dependency of the permeability 
in dependence on the relaxation time t and thus the fluid viscosity can be observed 
if standard bounce back boundary conditions are used. Using LB-MRT method the 
K-T correlation is significantly smaller than for LB-BGK l|20ll3T| . To keep focus in the 
accuracy comparison of the different setups, we use the same relaxation times for 
all presented calculations. The values are set to r/Ai = 0.857 and Tbuik/^^ = 0.84, 
because we found them to produce very accurate results for 3D Poiseuille flow in a 
previous contribution [32] . The first test geometry for measuring the permeability is a 
3D Poiseuille flow with different pipe lengths d/Ax. Fig.|4]shows the relative error of 
the computed permeabilities for chambers of different length and the three different 
setups presented. In case of p-BC and I-Ch the results are very accurate and the error 
stays below 2.5%. For Vp-S, however, there is a strong dependency of the results on 
the length of the chamber. If the sample is not periodic as it would be for experimen- 
tally relevant samples such as a discretization of a sandstone, it is not possible to avoid 
the chambers completely. In fact, the chambers are then absolutely needed providing 
a fluid reservoir before and after the sample, but as can be seen in Fig.|4] Vp-S, they 
produce a massive loss in accuracy. For this reason we do not take into account the 
Vp-S setup for further calculations. Our finding is of particular interest because this 
setup is the most popular one in the literature on permeability measurements using 
the LB method. 

On the other hand, for the results obtained by alternatives p-BC and I-Ch, one 
should take into account that the major error in permeability estimation of a real sam- 
ple with the LB method is due to the discretization of the porous space. For the 3D 
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Figure 4: Relative error of permeability of a 3D 
Poiseuille flow in a square pipe. Three different al- 
ternatives to drive the flow p-BC, Vp-S and l-Ch 
are tested together with different lengths of the 
chambers I and O. Using the very popular alterna- 
tive Vp-S, the results are highly dependent on the 
length of the chambers. Therefore, Vp-S is not 
taken into account for any further studies. 



Poiseuille flow in a square pipe this error is minimal since the pipe is aligned with the 
lattice. Therefore, a systematic error of the order of 2% is satisfactory. 

To estimate the effect of the chamber length on the permeability results when the 
Vp-S setup is being used, we simulate a 2D Poiseuille flow to compare the velocity 
field Ujj obtained inside the sample when chambers of two different lengths are used. 
To visualize the difference a relative error field e is calculated as 



e(x,A,B) 



u 



Xl 



(X,X; 



I/O 



A] 



u 



^1 



(X,X; 



I/O 



B) 



Ur 



(7.1) 



where A and B represents the two chamber lengths whose velocity field U;ci are be- 
ing compared. Fig. |5] displays two examples comparing the results of Xj/o = arid 
Xx/o = 2, and the results of Xj/o = 2 and Xj/o = 10. As can be clearly seen, the 
difference of permeability estimation when using different chamber lengths, can be ex- 
plained by the influence of the chamber geometry on the velocity fields u^j . Not only 
the velocity field entering the sample is dependent on the chamber length, but also the 
velocity field inside the sample, where the parabolic profile is already adopted. This 
issue is clearly shown in both examples of Fig. |5l where a constant relative difference 
in the velocity field u^j remains up to the fifth node inside of the sample. It is actually 
this difference, which responds to how the fluid enters the sample, that is responsi- 
ble for the dependency of the permeability on the length of the chamber. This clearly 
explains the results shown in Fig.H] Vp-S. 
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Figure 5: Two examples comparing the difference in the velocity field u^i obtained in a 2D Poiseuille flow 
d/Ax = 10 when the fluid is driven using Vp-S and two different chamber lengths Xj/q are used. The 
difference is calculated using the relative error field e{x,A,B) presented in Eq. i7.1\ . For both examples, the 
velocity field u^^j is only comparable in the zone that both domains have in common and which is defined 
by the shortest chamber. On the X2"a><is the first and last nodes inside the sample S represent the solid 
walls defining the pipe. Although the major difference resides at the inlet, a constant relative error remains 
inside the sample, which is responsible for the difference in the measured permeability. 

Fig.|6]displays the permeability calculated for a BCC cubic array using the LB-BGK 
and LB-MRT method and the alternatives p-BC and I-Ch with different chamber length. 
In contrast to the 3D Poiseuille flow, the measurements suffer indeed from an error due 
the domain discretization. As in Fig. H] for the 3D Poiseuille flow^, the dependency of 
the results using I-Ch on the length of the chambers is narrower than using p-BC for 
both LB methods. 

As it was stressed above, typically the major contribution to the error of perme- 
abilities estimated using the LB method is the discretization. This issue is clearly seen 
in Fig. [6] where the relative error in the permeability estimation increases if the radius 
of the sphere is decreased, which is equivalent to reducing the resolution. 

The LB-MRT method in general produces more accurate results than LB-BGK. The 
best results can be obtained when using LB-MRT together with I-Ch and Xj/o = 20. 
However, reducing the channel length does not have a strong influence. On the other 
hand, if we compare the results yield by using I-Ch with the ones obtained with p-BC, 
only minor differences can be obtained for identical chamber lengths. Following this 
we can assert that there is no significant difference in accuracy when using either I-Ch 
or p-BC, as it was also expected from the data presented in Fig. HI The parameter 
impacting on the accuracy is the chamber length. The LB-MRT together with p-BC 
without chambers achieves the same accuracy as LB-BGK with I-Ch, but in the latter 
case chambers are required. Since long chambers increase the computational domain, 
they are a burden in terms of computational effort. Thus, it is important to understand 
their influence in order to decide whether a lower accuracy can be accepted if the 
computational effort is substantially reduced. 

The efficiency of the LB-implementations, i.e., pore-matrix and pore-list as their 
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Figure 6: Relative error of permeability for a BCC cubic array of spheres. The alternatives p-BC, l-Ch 
are being tested together with different lengths of the chambers X and O and also using the LB-BGK and 
LB-MRT method. The LB-MRT together with p-BC without chambers achieves the same accuracy as using 
LB-BGK with l-Ch, but in the latter case chambers are required. 



memory demands explained above are directly related to the porosity of the sample. 
It is important to stress that the efficiency is independent of the setup used for the 
permeability estimation but it depends on the LB method used, i.e., BGK or MRT. 
In the case of using LB-BGK instead of LB-MRT the CPU time required decreases by 
15%-20% [i30il . 

The use of chambers and more precisely their lengths determine the porosity of 
the whole computational domain. For example, if in a cubic sample representing a 
realistic porous medium with porosity ^ = 10% we add chambers with lengths 5% of 
the sample length, we obtain a computational domain with porosity (p k, 18.2%. Even 
more dramatically is the case if the sample has a smaller porosity. If we add chambers 
with the same length explained above to a cubic sample with porosity (p = 5%, the 
total porosity would increase by a factor of ~173%. Since the efficiency of both LB- 
implementations scales linearly with the porosity, the use of chambers is expensive 
in terms of computational cost and should be avoided if possible. For low porosity 
values, as realistic samples have, the pore-list implementation is more efficient than 
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pore-matrix. On the other hand, for more open domains, as pipes for example, the 
pore-matrix data structure is more efficient since it does not require the connection 
matrix. A quantitative comparison of the computation efficiency is not easily possi- 
ble. Even though it is possible to calculate for both implementations the number of 
floating-point operations necessary to perform a single time step, a quantitative anal- 
ysis of the performance of two idependent codes is not as straightforward: the code 
performance strongly depends on the actual implementation details, memory access 
patterns and in particular also on the cache performance of the used architecture. 



8 Conclusions 

We presented different alternative setups to estimate the permeability of a porous 
medium, which utilize an injection channel (I-Ch), pressure boundary conditions (p-BC), 
or a force density applied in the sample volume (Vp-S). The setups differ on how the 
fluid is driven and how the pressure gradient is being estimated. Further, the accu- 
racy of the LB-BGK and LB-MRT method together with different setups with variable 
chamber length, has been investigated. Another point that has been stressed is the 
efficiency of the LB-implementations comparing the standard implementation using a 
pore-matrix data structure with a pore-list code to represent the geometry. 

Our findings show that the use of Vp-S as an alternative for permeability estima- 
tion should not be taken into account due to the need of chambers as fluid reservoir 
together with the strong dependency of the results on the length of these chambers. 
This is an important result since driving the fluid by using a body force applied in the 
sample volume is probably the most popular approach in the literature on permeabil- 
ity calculations using the LB method. 

Further, the pore-list implementation allows to reduce the computational effort 
required to simulate fluid flows in porous media substantially if the porosity of the 
porous medium is small. 

Finally, taking into account the extra computational effort added because of the 
number of required lattice nodes when chambers are used, it is highly recommended 
to completely eliminate them by using p-BC. A resulting compromise in the accuracy 
can be resolved by using the LB-MRT method. For this reason in the case of facing 
the permeability estimation of a realistic media, pore-list, p-BC, and LB-MRT is the 
combination suggested by the authors taking into account the computational domain 
size, accuracy and required computer time. 
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